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ABSTRACT 

We assess the probable redshift {zj-^ ~ 7) for full reionization of the inter- 
galactic medium (IGM) using a prescription for the co-moving star-formation- 
rate (SFR) density (psfr) required to maintain photoionization against recom- 
bination. Our newly developed on-line reionization simulator allows users to 
assess the required SFR and ionization histories, using a variety of assumptions 
for galactic and stellar populations, IGM clumping factor and temperature, and 
LyC escape fraction. The decline in high-redshift galaxy candidates and Lya 
emitters at 2 = 6 — 8 suggests a rising neutral fraction, with reionization at 
z^7 increasingly difficult owing to increased recombination rates and constraints 
from the ionizing background and LyC mean free path. The required rate is 
Psfr ~ (0.018 Mq yr'^ Mpc-^)[{1 + z)/8f{CH/S){0.2/ f,,,)T^''-'''' scaled to fidu- 
cial values of clumping factor Ch = 3, escape fraction /esc = 0.2, electron tem- 
perature Te = 10^ K, and low-metallicity initial mass functions (IMF) and stellar 
atmospheres. Our hydrodynamical + N-body simulations find a mean clumping 
factor Ch ~ (2-9)[(l + z)/Q]~^'^ in the photoionized, photoheated filaments at 
z = 5 — 9. The critical SFR could be reduced by increasing the minimum stellar 
mass, invoking a top-heavy IMF, or systematically increasing /esc at high z. The 
CMB optical depth, Tg = 0.088 ± 0.015, could be explained by full reionization, 
producing = 0.050 back to Zj-^i ~ 7, augmented by Are ~ 0.01 — 0.04 in a par- 
tially ionized IGM at 2 > 7. In this scenario, the strongest 21-cm signal should 
occur at redshifted frequencies 124 — 167 MHz owing to IGM heating over an 
interval Az ^ 3 from z ^ 7.5 — 10.5. 
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Subject headings: Galaxies: high-redshift — intergalactic medium — reionization 

1. INTRODUCTION 

Our current hypothesis is that the galaxies and black holes observed today originated 
over 13 Gyr ago, growing from seeds of primordial density perturbations. One can test this 
hypothesis by studying the star formation rate (SFR) history, from the Epoch of Reionization 
(EoR) at redshifts 6 < z < 12, through the peak era at z ^ 2 — 3 (Hopkins & Beacom 2006) 
down to the present epoch. The feedback of ionizing radiation, kinetic energy, and heavy 
elements leaves imprints on early stars, supernovae, and galaxies, providing a "fossil record" 
that can be detected through abundances in Galactic halo stars and the intergalactic medium 
(IGM) and in the distributions of mass, metallicity, and luminosity of galaxies. 

Determining when and how the universe was reionized by these early sources have been 
important questions for decades (Gunn & Peterson 1965; Sunyaev 1977; Robertson et al. 
2010). It has been suggested that IGM reionization was complete by z ^ 6.5 (Fan et al. 2001; 
Gnedin & Fan 2006; Fan et al. 2006; Hu & Cowie 2006), based on strong Lya absorption from 
neutral hydrogen along lines of sight to QSOs at z > Q. Becker et al. (2007) and Songaila 
(2004) used transmission of the Lya (and Ly/3) forest out to z = 5.8 and z = 6.3, respectively, 
to suggest a smoothly decreasing ionization rate toward higher redshifts. Recent surveys of 
high-redshift galaxies and Lya emitters (Bouwens et al. 2011a; Ouchi et al. 2010; Kashikawa 
et al. 2011; Ono et al. 2012; Schenker et al. 2012) infer an increasing IGM neutral fraction 
from the declining populations between z = Q — 8. Further evidence for an increasing neutral 
fraction comes from the decreasing sizes of ionized "near zones" associated with quasars 
between z = 5.7 and z = 6.4 (Carilli et al. 2010) and from the Lya damping wing in the 
transmission profile toward the newly discovered quasar at z = 7.085 (Mortlock et al. 2011). 
These studies all suggest that the IGM is becoming increasingly neutral between z = 6 — 7, 
marking the end of cosmic reionization when ionized regions overlap and percolate. Whether 
the epoch of full reionization occurs at Zj-^ ~ 7 is still not ascertained. 

A contrasting estimate of the EoR comes from the measured optical depth, Tg = 
0.088±0.015, to electron scattering of the cosmic microwave background (CMB) in WMAP-7 
observations (Larson et al. 2011; Komatsu et al. 2011). The error bars come from the central 
68% in the marginalized cumulative distribution. Using additional cosmological parameter 
constraints, they infer single-epoch reionization at Zj-d = 10.6 ± 1.2. Although such a high 
redshift could be explained with ACDM simulations and modeled SFR histories (Choudhury 
& Ferrara 2005; Trac & Gen 2007; ShuU & Venkatesan 2008), the GMB observations are at 
variance with optical surveys that suggest late reionization, unless reionization is a process 
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that extends to higher redshifts. 

For a fully ionized IGM, including both H and He, the optical depth r^. = 0.050 for 
^rei = 7.0 (see Section 2.1). In the CMB analysis (Larson et al. 2011), marginalization of 
Tf. with other cosmological parameters allows the possibility of lower optical depth, with a 
reionization epoch as low as 2rci ~ 7 — 8 at 95% CL. One can also invoke a partially ionized 
IGM at Zj-ci > 7, as discussed by many groups (Cen 2003; Venkatesan et al. 2003; Ricotti & 
Ostriker 2004; Benson et al. 2006; Shull & Venkatesan 2008). However, even with the recent 
progress in finding high-^ galaxies, we still do not know whether galaxies are the sole agents of 
reionization. Current observations of high-z galaxies leave open several ionization scenarios, 
some involving simple hydrogen reionization at z «i 10 and He II reionization at z ^ 3, 
and others with more complex ionization histories (Bolton & Haehnelt 2007; Venkatesan, 
Tumlinson & Shull 2003; Cen 2003) that depend on SFRs at ^ = 7 - 20. Shull & Venkatesan 
(2008) demonstrated how the CMB optical depth constrains the SFR and IGM metallicity 
history at 2 > 7, and Trenti & Shull (2010) quantified the metallicity-driven transition from 
Population HI (metal-free stars) to Population II (stars formed from metal-enriched gas). 

Several recent observations provide valuable constraints on the luminosity function of 
high- 2; galaxies. The number density of galaxies appears to drop rapidly at z > 7 (Bouwens 
et al. 2009, 2010a,b, 2011a,b,c). With a comoving SFR density p*<0.01 Mq yr-^ Mpc"^ 
at z ^ 7 (Gonzalez et al. 2010; Bouwens et al. 2011a), the observed galaxies do not pro- 
duce enough ionizing photons in the Lyman continuum (LyC) to maintain a photoionized 
IGM against recombinations. However, the luminosity function is steep, and the total LyC 
budget requires extrapolation to low-luminosity galaxies (Trenti et al. 2010; Bouwens et al. 
2011b). Moreover, the conversion from SFR to LyC production rate relies on insecure cali- 
brations from theoretical models and comparison with high-mass, low-metallicity stars. We 
revisit the calculation of LyC photon production and assess the high- 2; galaxy contribution 
to reionization. We also analyze several factors, such as the photon escape fraction (/esc); 
IGM clumping factor (Ch), and electron temperature (Te), which enter the calculation of 
the "critical star formation rate" (pcrit) necessary to maintain a photoionized IGM. 

In Section 2, we calculate pcrit {Mq jt~^ Mpc~^) in a filamentary IGM, equating the 
production rate of Lyman continuum (LyC) photons with the hydrogen recombination rate. 
The photoionization rate depends on the mass function of stellar populations, their evolu- 
tionary tracks and stellar atmospheres, and the escape fraction, /esc of LyC photons away 
from their galactic sources. The recombination rate depends on the density and tempera- 
ture of the IGM, properties we explore with cosmological simulations. In Section 3, we give 
our results for the critical SFR at z > 6 and present our new SFR simulator, a user-friendly 
interface for calculating pcrit(^) and Te{z). In Section 4, we discuss the implications for the 
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hydrogen EoR. Consistency between high-redshift galaxies and CMB optical depth appears 
to require Zt-^i ~ 7 and a partially ionized IGM at 2; > 7. The peak signal from redshifted 
21-cm emission would likely occur during the heating period between z = 7.7 — 8.8 (145-163 
MHz) when the hydrogen neutral fraction xhi ~ 0.5 (Pritchard et al. 2010; Lidz et al. 2008). 



REIONIZATION WITH CLUMPING AND PHOTON ESCAPE 



2.1. Reionization and Critical Star Formation Rate 



We denote by Psfr (^0 yr~^ Mpc^^) the global star formation rate per co-moving 
volume. Using a simple argument (Madau et al. 1999), balancing photoionization with 
radiative recombination, we estimate the critical SFR to maintain IGM photoionization at 
z > 7, assuming that the LyC photons are produced by populations of massive (OB-type) 
stars. Because the mass in collapsed objects (clusters, groups, galaxies) is still small at high 
redshift, the IGM contains most of the cosmological baryons, at mean density 

Pb = QbPcril + zf = (4.24 X 10-31 g cm-3)(l + zf . (1) 

For a Hubble constant denoted Hq = (100 km s"^ Mpc"'^) h, we adopt the WMAP-7 (plus 
BAG + Ho) parameters, Qbh^ = 0.02255 ± 0.00054 and Q^h^ = 0.1352 ± 0.0036 (Komatsu 
et al. 2011) relative to a critical density pcr = 1.8785 x 10~^^ h"^ g cm"^. From the corre- 
sponding helium mass fraction Y = 0.2477 ±0.0029 (Peimbert et al. 2007), we adopt a mean 
hydrogen number density, 

Ptil - Y) 



rriH 



;i.905 X 10"' cm-^)(l + 2) 



(2) 



In a fully ionized IGM, the CMB optical depth back to be written as the integral 

of UgaTdi, the electron density times the Thomson cross section along proper length. 



TAZr 



Hearil + z)~^ [c/H{z)] dz 



(3) 



for a standard ACDM cosmology {Qm + Qa = 1) with H{z) = i/o[^m(l + + QaY^"^. This 
integral can be done analytically (Shull & Venkatesan 2008), 



''"e(^rei) 



CCTt 



2Qb 



p,,{l-Y){l + y) 



In the high-redshift limit, when Q 



rriH 



[{fi^(l + 2rei)' + ^^A}'/'-l] . (4) 



[l + z)^^ this expression simplifies to 



TAZr, 



2Qb 



p,,{l-Y){l+y) 



[l + z, 



.^3/2 



(0.0521) 



:i + z. 



3/2 



(5) 
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independent of h to lowest order. The helium and electron densities are written tt-hg = ynn 
and TT-e = nH{l + y) for singly ionized helium, where y = rinc/nu = (F/4)/(l — F) ^ 0.0823 by 
number. To these formulae, we add Are ~ 0.002, from electrons donated by He III reionized 
at z < 3 (ShuU et al. 2010). Helium therefore contributes ~8% to t^, and a fully ionized 
IGM produces Tg = 0.050, 0.060, and 0.070 back to redshifts z^ei = 7, 8, and 9, respectively. 

A comoving volume of 1 Mpc^ contains A'^^ = 5.6 x 10^^ hydrogen atoms. Our simple ion- 
ization criterion requires a SFR density that produces a number of LyC photons equal to Nh 
over a hydrogen recombination time, tree = ['^eCfc/f The hydrogen Case-B recombina- 

tion rate coefficient (Osterbrock & Ferland 2006) is a5f^(T) ^ (2.59 x lO'^^ ^^^3 ^-^-^t^o-^^^ ^ 
scaled to an IGM temperature T = (10^ K)T4. For typical IGM ionization histories and pho- 
toelectric heating rates, numerical simulations predict that diffuse photoionized filaments of 
hydrogen have temperatures ranging from 5000 K to 20,000 K (Dave et al. 2001; Smith et al. 
2011). These are consistent with temperatures inferred from observations of the Lya forest 
at 2 < 5 (Becker et al. 2011). 

Owing to gravitational instability, a realistic IGM is inhomogenous and filamentary. 
Semi-analytical models of the reionization of the universe often adopt a "clumping factor" , 
Ch = (nl) / (ne)"^, to account for inhomogeneity in estimates of the enhanced recombination 
rate in denser IGM filaments. The clumping factor therefore plays an important role in 
computing the critical SFR density needed to maintain the reionization of the universe. The 
clumping factor is also used in numerical simulations to implement "sub-grid physics", in 
which changes in the density field occur on scales below the resolution of the simulation 
and are also approximated by the factor Ch (Gnedin & Ostriker 1997; Madau et al. 1999; 
Miralda-Escude et al. 2000; Miralda-Escude 2003; Kohler et al. 2007). 

We correct the recombination time for density variations scaled to a fiducial Cj^ ~ 3, 
found in the simulations described below. At z ~ 7, the IGM filaments have electron density 
rif. ^ (10~^ cm~^)[(l + z)/8]^Ch, and the characteristic times for hydrogen recombination 
and Hubble expansion are, 

(6) 

tn ^ [HoQ'd'il + zf/'r'^il.lSGyi) 

In our calculations, we express the reionization criterion as A'i^(Mpc~^) = Pcrit^rocQLyc/osc, 
where pcrit is the critical SFR density {Mq yr~^ Mpc""^) and QljC is the conversion factor 
from psFR to the LyC production rate (see Section 2.2). We define /esc as the fraction of 
LyC photons that escape from their galactic sources into the IGM (Dove & Shull 1994). 



tree ~ (386 Myr)(3/CH)T° 



-o/ z 



(7) 
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Recent statistical estimates (Nestor et al. 2011) suggest that /esc ~ 0.1 for an ensemble of 
26 Lyman-break galaxies and 130 Lya emitters at z ~ 3.09 ± 0.03, and it could be higher 
for the lower-mass galaxies that likely dominate the escaping LyC at z > 6 (Fernandez & 
ShuU 2011). The LyC production efficiency, QhyC: is expressed in units 10^^ photons per 
Mq of star formation, since typical massive stars emit (1 — 10) x 10^^ LyC photons over 
their lifetime. To evaluate QhyC, we convert the SFR (by mass) into numbers of OB-stars 
and compute the total number of ionizing photons produced by a star of mass m over its 
lifetime. We then integrate over an IMF, \l/(m) = Km""', with a range m^[^ < m < m^^^. 
The standard mass range is 0.1 Mq to 100 M©, but changes in the mass range and IMF 
slope will affect the LyC production substantially. For example, differences in the IMF have 
been associated with higher Jeans masses in low-met allicity gas in the high-redshift IGM 
(Abel et al. 2002; Bromm & Loeb 2003), and Tumlinson (2007) and Smith et al. (2009) 
noted the potential influence of CMB temperature on modes of high-redshift star formation. 
Consequently, most cosmological simulations or calculations include a metallicity-induced 
IMF transition (Trenti & Shull 2010) between low-met allicity Population III star formation 
and metal-enhanced Population II. 



2.2. Ionizing photon luminosities 

The number of LyC photons produced per total mass in star formation depends on the 
IMF of the stellar population and is given by the conversion factor QLyC- We calculate this 
conversion by integrating the total number of LyC photons produced over the entire mass 
in star formation. 

Nt n f^""''' "^(171)0(171) dm 
PsFRtrec L , "i! [m)m dm 

Here, Q{m) is the lifetime-integrated number of LyC photons as a function of mass calculated 
from stellar atmosphere models and evolutionary tracks. The IMF, \l/(m) = Km'"', is 
integrated over the mass range m^[^ < m < mmax? where m is expressed in solar units. 
In the SFR simulator discussed in Section 3.3, the user can choose between a normal and 
broken IMF power law. The lower integration limit in the numerator, moB, is the mass at 
which stars no longer produce significant amounts of LyC photons. 

In our calculator, we use and compare two models that calculate Q{m). First, using 
stellar atmospheres and evolutionary tracks (R. S. Sutherland & J. M. Shull, unpublished), 
we find that, over its main-sequence and post- main- sequence lifetime, an OB star of mass 
m produces a total number iVLyc = Q63("^) x 10^^ of ionizing photons, where Qes ~ 1 — 10 
over the mass range m = 30 — 100 and for metallicities Z = (0.02 — 2.0)Zq. We have fitted 
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our results to the form Qesim) ^ Am — B, for m > moB = B/A (the mass moB defines 
the effective lower limit for stars that produce significant numbers of LyC photons). For 
metallicities in the range 0.002 < Z < 0.04 (where Z = Zq = 0.02 is the solar metallicity) 
and for masses 15 < m < 100, the fitted coefficients are B = 1.578 and A{Z) = 0.0950 — 
1.059Z. Thus, for Z = Zq ^ 0.02, we have A = 0.0738 and moB = 21.4 Mq, while for 
Z = O.IZq = 0.002 we have A = 0.0929 and moB = 17 M©. Inserting these values for Q{m) 
into equation (8), we integrate these LyC photon yields over the IMF, to derive the conversion 
coefficient, QhyC: from total mass in star formation to total number of LyC photons produced 
(in units of 10^'^ photons per Mq ). 

The integrals in eq. (8) can be done analytically as functions of the IMF parameters 
and metalhcity. For a Salpeter IMF (a = 2.35, 0.1 < m < 100) we find: QtyC = 0.00236 (for 
Z = Zq), QLyC = 0.00366 (for Z = 0.2Zq), and QljC = 0.00384 (for Z = O.IZq). Adopting 
a typical (low-Z) value Qhyc = 0.004 (4 x 10^° photons per Mq), we can solve for the critical 
SFR for reionization, scaled to clumping factor Ch = 3, escape fraction /esc = 0.2, and gas 
temperature T4 = 1: 



Pcrit = (0.018 Mq yr-1 Mpc" 



+ 



L/esc/0.2j 



0.004 



-0.845 



(9) 



Our chosen value of Cji^ ~ 3 is consistent with recent downward revisions (Pawlik et al. 2009) 
and with our numerical simulations discussed in Section 2.3. Escape fractions /esc ~ 0.1 — 0.2 
have been inferred from observations at 2; ~ 3 (Shapley et al. 2006; Nestor et al. 2011) 
and theoretical expectations (Fernandez & Shull 2011). Equation (9) agrees with prior 
estimates (Madau et al. 1999) when adjusted for our new scaling factors, in particular the 
ratio {Ch If CSC = 15). Their earlier paper assumed Ch = 30, /esc = 1, ^bh'^ = 0.020, 
^LyC = 0.005, and z = Our fiducial redshift has increased from z = 5 to 2; = 7, appropriate 
for the new discoveries of high- redshift galaxy candidates for reionization. Expressed with the 
same coefficients in eq. (9), their coefficient would be nearly the same, 0.020 Mq yr~^ Mpc~^ 
at z = 7. One of the improvements in our formulation is to better identify the dependences 
on the physical parameters {Ch, /esc, T^) and the SFR-to-LyC conversion factor (QLyc), 
which can change with different IMFs and atmospheres. 

A related calculation is the production rate of ionizing (LyC) photons per unit volume, 
needed to balance hydrogen recombinations. With the same assumptions as above, this is 



LyC 



dt 



n, 



nnii ajf ^ = (4.6 x 10^° s"^ Mpc" 



crit 



-0.845 



For standard high-mass stars (07 V, solar metallicity) each with LyC photon luminosity 
10^^ s~^, this rate corresponds to an 0-star space density uq ~ 50 Mpc"^ at 2; = 7. Schaerer 
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(2002, 2003) also computed models of Population III and low-metallicity stars based on non- 
LTE model atmospheres and new stellar evolution tracks and evolutionary synthesis models. 
The lifetime total number of LyC photons produced per star of mass m is Q{m) = Q{H) x t^. 
For stars of mass parameter x = log{m/MQ), the number of ionizing photons, Q, emitted 
per second per star, is given by: 



logio [Qh/ s ^] 



43.61 +4.90X- 0.83x2 ^ 

39.29 + 8.55X Z 

27.80 + 30.68X - 14.80x2 + 2.50x3 Z 

27.89 + 27.75X- 11.87x2 + 1.73x3 Z 



0, 
0, 

O.O2Z0, 



9 - 5OOM0, 
5 - 9Mq, 
7 - ISOMq, 
7 - I2OM0 



and the star's lifetime is given by: 



r 9.785 - 3. 759x + 1.413x2 - 0.186x3 Z = 0, 
logio [i*] = { 9.59 - 2.79x + 0.63x2 ^ ^ 0.02Zo, (12) 

[ 9.986 - 3.497X + 0.894x2 Z = Zq 

Sample values of Qljc and pcrit for different IMFs are shown in Table 1, for both model 
atmospheres. One can obtain different values of QljC for power-law IMFs by varying their 
high-mass slope /3 and the minimum and maximum masses. Increases in Qhyc translate into 
decreases in critical SFR. For example, at 0.1 Z©, if the minimum mass is raised to 1 M© 
with a Salpeter slope (/3 = 2.35), the LyC production factor rises to QhyC = 0.0098, some 
2.5 times higher than the equivalent model, QhyC = 0.00384 for m^[^ = O.IMq. If the IMF 
is flatter, /3 = 2 instead of 2.35, with minimum mass fixed at m^[^ = O.IMq, QljC = 0.0177 
(4.6 times higher). At the upper end of the IMF, if one increases mmax from 100 to 200 Mq, 
keeping (3 = 2.35 and mmin = 0.1, one finds QljC = 0.0054 (1.4 times higher). All of these 
variations can be explored with our reionization/SFR simulator, described in Section 3.3. 

The critical SFR of 0.018 Mq yr~^ Mpc"^ can be understood from simple arithmetic. 
Over the 386 Myr recombination time (eq. [6]) with Ch = 3, /esc = 0.2, and z = 7, approxi- 
mately 20 million stars are formed per comoving Mpc^ in a Salpeter IMF (0.1 — 100 Mq) with 
mean stellar mass O.352M0. Of these stars, a small fraction, /qb ~ 7x 10~^, are massive OB 
stars (m > 2OM0). At an average of 2 x 10^^ LyQ photons and escape fraction f^sc = 0.2, 
these ~14,000 OB stars will produce a net (escaping) ~ 5 x 10^^ LyC photons over their 
lifetime. These photons are sufficient to ionize Nh = 5.6 x 10^^ hydrogen atoms Mpc"^. 
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2.3. Clumping Factor from Cosmological Simulations 

The simulations used in tliis study were performed with Enzo, an Eulerian adaptive 
mesh-refinement (AMR), hydrodynamical + N-body code (Bryan & Norman 1997; O'Shea 
et al. 2004, 2005). Smith et al. (2011) enhanced this code by adding new modules for star 
formation, primordial chemistry, and cooling rates consistent with ionizing radiation, metal 
transport, and feedback. The ionizing background is spatially constant and optically thin, 
but variable in redshift. At this stage, we have not implemented radiative transfer or spectral 
filtering by the IGM. 

For the clumping factor calculation, we ran a simulation on a 50h~^ Mpc static grid 
("unigrid") cube with 1024^ cells, denoted as run 50_1024_2 in Table 1 of Smith et al. (2011). 
The radiative heating from the ionizing background plays an important role in determining 
the properties of the filamentary structure. As a filament is ionized and heated, its density 
drops and its temperature rises; both effects reduce the recombination rate. To study these 
effects, we ran four moderate-resolution [50h~^ Mpc unigrid cube with 512^ cells) simulations 
with varying ionizing backgrounds, summarized in Table 2. After initial submission of this 
paper, we ran a 1536'^ simulation, to check convergence and assess the "cosmic variance" 
among eight sub-volumes of the 1024^ and 1536^ simulations. The standard UV background 
was taken from Haardt & Madau (2001), although we also explore new computations of 
high- 2; SFRs by Trenti et al. (2010) and Haardt & Madau (2012). These four simulations 
were: (1) no photoionizing background; (2) UV background ramped up from z = 7 to z = 6 
(run 50_512_2 from Smith et al. 2011); (3) UV background ramped up from z = 9 to z = 8; 
and (4) UV background ramped up from z = 9 to z = 8, but twice as strong as in (3). 
Post-processing of the simulations was performed using the data analysis and visualization 
package, yt§, documented by Turk et al. (2011). 

For regions of ionized hydrogen of density nnii, we calculate the clumping factor, Ch- 
using two different methods. The first calculation, which has been used in some earlier 
studies, uses density weighting. In this "density field" (DF) method we define 

Cdf = j ^"A , (13) 

and average the density fields (quantities Xi, denoting either nnii or n^jj) where parameter 
averages are computed by summing over grid cells (j), subject to various "cuts" on the IGM 



^http: / /yt. enzotools.org/ 



overdensity and gas temperature and weighted by factors, wj, 



(14) 



j j 



In the second calculation, we compare the local recombination rate to the global average 
recombination rate, averaged over density and temperature in cells. 



where again (x) denotes a weighted average and a^j (T) is the case-B radiative recombination 
rate coefficient for hydrogen, as tabulated by Osterbrock & Ferland (2006). We refer to this 
as the "recombination rate" (RR) method. 

Because the purpose of the clumping factor is to correct for an enhanced recombination 
rate, we believe Crr to be a more appropriate representation. The recombinations that 
are important to removing ionizing photons occur in the filamentary structure of the IGM, 
and the clumping factor should only be calculated in these structures. To assess the critical 
SFR necessary to maintain an ionized medium, we focus on grid cells that are significantly 
ionized. Because a negligible amount of recombination occurs in cells containing mostly 
neutral gas, these cells should not contribute to the clumping factor. If we do not exclude 
neutral gas, the clumping factor is extremely high {Ch ~ 100) before the ionizing background 
is turned on. We find large density gradients between the neutral gas (uniformly distributed 
over the simulation) and the ionized gas. Because the clumping factor is a measure of 
the inhomogeneity of the medium, a large density gradient yields a large value of Ch- If 
low-density voids are included in the calculation, the larger density gradient leads to an 
overestimate of the clumping. By setting both upper and lower density thresholds, we can 
exclude collapsed halos and low- density voids from our calculations. 

Previous studies (Miralda-Escude 2000; Miralda-Escude et al. 2003; Pawlik et al. 2009) 
addressed these issues by setting a density threshold that excludes collapsed halos, but they 
did not set a lower limit to exclude voids from their calculations. To explore these effects in a 
ffiamentary IGM, we make various cuts of our data in baryon overdensity (A;, = pb/ pb) and in 
temperature, metallicity, and hydrogen ionization fraction (x = nYui/nH). In our standard 
formulation, we include only those cells that meet the following criteria: 1 < < 100, 
300 K < T < 10^ K, Z < lO"'' Zq, and x > 0.05. We believe this "data cut" adequately 
represents unenriched IGM lying on an adiabat (see Figure 19 of Smith et al. 2011). We 
also explore the clumping factor with no lower density threshold (total range A;, < 100). 
Our results, presented in Section 3.1, show a small increase in Ch{z) from the wider range 



K)(nHii) («i?^(T)) ' 




(15) 
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in densities by including low-density cells with Af, < 1. However, these low-density voids do 
not contribute substantially to the recombination rate. 

In summary, our prescriptions for calculating the clumping factor yield a more physi- 
cal representation of the enhanced recombination rate that is an important component to 
many reionization models. By not assuming a fully ionized medium and by specifically fol- 
lowing nmi, we are able to exclude denser neutral gas that does not contribute appreciably 
to the recombination rate. We make simple assumptions on the reionization process and 
reionization history (Section 3.2), turning on the ionizing radiation field at redshifts 2; = 7 or 
z = 9 and following the thermal history arising from photoelectric heating, radiative cooling, 
and pressure smoothing (Pawlik et al. 2009). The metal-line and molecular cooling, metal 
transport, and feedback included in our simulations also allow us to accurately represent 
the thermodynamics of the gas, which has been shown to have a significant effect on the 
evolution of the clumping factor. Future models will include discrete sources and radiative 
transfer, accounting for temperature increases arising from photo-heating with spectral hard- 
ening (Abel & Haehnelt 1999). Because our current simulations employ a spatially constant 
ionizing radiation field, we anticipate carrying out these more realistic situations. 



3. RESULTS FROM SIMULATIONS 

3.1. Computing the Clumping Factor 

Early studies of IGM clumping adopted high values, Ch > 40 at 2; < 5 (Gnedin & 
Ostriker 1997). As discussed earlier, we believe these values are too high for the ionized 
IGM filaments, which have expanded as a result of the heat deposited by LyC photons. The 
differences between observations and inferred critical SFR densities can largely be attributed 
to this high clumping factor (Sawicki & Thompson 2006; Bouwens et al. 2007). More recent 
studies have trended towards less clumping. Raicevic & Theuns (2011) argue that using a 
global clumping factor overestimates the recombination rate, and that local values should be 
used instead. In this study, we calculate the clumping factor for a series of high-resolution 
cosmological simulations for ionized hydrogen and helium that explore how the photoheating 
of an ionizing background affects the IGM thermodynamics (density and temperature). We 
impose criteria in overdensity, temperature, metallicity, and ionization fraction to constrain 
our calculations to IGM filaments where recombinations and ionization fronts are most im- 
portant. The clumping in simulation 50_1024_7 is calculated as a function of redshift from 
z = 15 to z = foT both Cdf and Crr, weighted by overdensity and by volume. Note 
that the unigrid simulations automatically yield volume weighting, when averaged over cells. 
For the remainder of our simulations, we find a power-law overdensity distribution, with an 
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average form /{Ah) = 10^ AjJ"^'^. Weighting these cells by overdensity gives undue emphasis 
on small numbers of high-density cells. This method is also biased, since it does not calculate 
the clumping factor where most of the recombinations are occurring. 

Figure 1 compares the DF and RR methods and shows the difference between the dif- 
ferent weights used in our calculations. In both methods, weighting by overdensity results 
in higher clumping factors. The recombination-rate method (eq. [15]) results in a lower 
clumping factor because it accounts for both density and temperature effects from the ion- 
izing background on the clumping averages. Photoelectric heating during photoionization 
causes the filaments to expand, lowering the density. The increased electron temperature 
also reduces the radiative recombination rate coefficient. The DF method (eq. [13]) includes 
no dependence on the recombination rate coefficient, and therefore only counts the effect of 
photoheating on the filament density. See Section 3.2 for further discussion. 

In our current simulations, we turn on a spatially constant UV background at 2 = 7 or 
z = 9. The effects of different "density cuts" in the summation over cells (see Section 2.3) 
are shown in Figure 2. We find little difference at 2; > 9 (turn-on of ionizing radiation), and 
see a small increase in Ch at 2 < 9 when we include the lower density cells with < 1. 
In future, more refined simulations with radiative transfer, we expect photoionization to be 
initiated in high-density regions, where the stars are located. Once the UV background is 
turned on, lower-density regions are easily photoionized. The photoelectric heating of the UV 
background causes the filaments to expand and become diffuse, resulting in a lower clumping 
factor. The elevated temperature also reduces the recombination rate coefficient, which in 
turn lowers Ch when calculated by the RR method. When calculating a mean recombination 
rate, only the recombining ionized gas is relevant. By 2; ~ 6, the IGM is nearly completely 
ionized and only small regions of H I remain. As noted earlier, in more refined simulations, 
the remaining H I would be largely self-shielded from the UV background and likely to reside 
in clumps of high density yet unreached by ionization fronts. 

Figure 3 explores the convergence and cosmic variance among different simulations. In 
two panels, we compare results from our 512^ simulations with larger simulations with 1024^ 
and 1536^ cells. In each panel, we show the average clumping factors, together with those 
in eight sub-volumes. The average values of Ch agree in the 1024^ and 1536^ simulations, 
and the small (±10%) variance among these sub- volumes suggests that the 1024^ and 1536^ 
simulations are converged. Over the redshift range 5 < 2; < 9 in the 1536^ simulation, the 
clumping factor is well fitted by a power-law, 

-1.1 

(16) 

representing a slow rise in clumping to lower redshift, after the turn-on of ionizing radiation. 



Ch{z) = (2.9) 
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3.2. Ionizing Background Study 

To study the effect of radiative lieating from tfie ionizing background on filamentary 
structure, we ran a suite of mo derate- resolution simulations witli different ionizing back- 
grounds (see Table 2 for details). The redshift at which the UV background is turned on 
affects the clumping factor, causing it to drop and then recover at lower redshift. Pawlik et al. 
(2009) attribute this effect to Jeans filtering, where the photo-heating raises the cosmological 
Jeans mass, preventing further accretion onto low-mass halos and smoothing out small-scale 
density fluctuations. The photoheating also heats the filaments we are focusing on, which 
lowers the hydrogen recombination rate coefficient, oc and causes the filaments 

to expand; both of these effects reduce the clumping factor. Without a UV background, the 
clumping continues to rise as the filaments gravitationally collapse. Pawlik et al. (2009) also 
claim that the clumping factor at z = 6 is insensitive to when the background is turned on, 
as long as it is turned on at 2; > 9. When an ionizing background is introduced, photoheating 
acts as a positive feedback to reionization by lowering the clumping factor and making it 
easier to stay ionized. This same photoheating mechanism also suppresses star formation in 
low-mass halos, which in turn lowers the ionizing photon production rate by star-forming 
galaxies and acts as a negative feedback to reionization (Pawlik et al. 2009). These thermo- 
dynamic processes affect clumping and structure and emphasize the importance of carefully 
modeling the strength of feedback processes and their effects on Jeans mass. 

We have explored what happens when the ionizing radiation turned on earlier, especially 
during the interval 2 = 6 — 8 marking the transition from a neutral to fully ionized IGM. In 
a series of 512^ simulations, we explore the influence of turn-on of photoionizing radiation 
between redshifts z = 7 and 2 = 9. Figure 4 shows Ch, computed for the moderate- 
resolution simulations via the density-field method and weighted by volume. We do not 
plot simulation 50_512_9_2, since it is identical to simulation 50_512_9. As in the results of 
Pawlik et al. (2009), we find that photoheating from the ionizing background results in a 
decrease in the clumping factor. The clumping factors fall along two tracks: a higher track 
at redshifts above the turn-on of the UV background, and a lower one after turn-on. By 
2; = 5, we find nearly identical clumping factors for both backgrounds. After a substantial 
recovery time, the redshift when the background is turned on is not important. Before this 
recovery, the clumping factor of the earlier background is lower by a factor of ~2, resulting 
in earlier reionization. 

One can compare the strengths of feedback to reionization, where "positive feedback" 
lowers the clumping factor and "negative feedback" suppresses star formation. At 2 = 5 
the total stellar mass (SFR density) of simulation 50_512_7 is 1.19 (1.13) times lower than 
that of simulation 50_512_0, while the clumping factor is 1.64 times lower. For simulation 
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50_512_9, the total stellar mass (SFR density) and clumping factor are 1.52 (1.49) and 1.66, 
respectively, times lower than those of simulation 50_512_0. This suggests that the positive 
feedback introduced by the background is greater than the negative feedback. However, the 
stellar mass (SFR density) does not recover in the same manner as the clumping factor. 
Once photoheating suppresses the formation of small-mass halos, the Hubble flow takes over 
and prevents them from collapsing and forming stars. Therefore, the redshift at which the 
background is turned on determines whether the positive or negative feedback dominates 
and whether the ionizing background will cause reionization to be accelerated or delayed. 



3.3. Reionization SFR Simulator 

In connection with this project, we have developed a user interface for calculating the 
critical SFR density (pcrit) needed to maintain the IGM ionization at a given redshift. The 
software computes the effects of variations in the stellar IMF (slope, mass-range) and model 
atmospheres, and the redshift evolution of metallicity and gas thermodynamics (density, 
temperature, coupling to the CMB). The clumping factor and LyC escape fraction are free 
parameters in the calculator. Our simulations find ranges of Ch ~ 1 — 10 depending on 
redshift, overdensity, and thermal phase of the (photoionized or shock-heated) IGM. For the 
new 1536^ simuation, the global mean clumping factor is {Ch) ~ 3, with a power-law fit 
Ch{z) = (2.9) [(1 + z)/6]-^-^ for redshifts between 5 < ;z < 9. 

This calculator is a useful tool for determining the population of galaxies responsible 
for reionization. The critical SFR per co- moving volume (Eq. [9]) is obtained by balancing 
the production rate of LyC photons with the number of hydrogen recombinations. Here, Ch 
is the clumping factor of ionized hydrogen, /esc is the escape fraction of LyC photons from 
their host galaxies, and T4 is the temperature scaled to 10^ K. The conversion factor, QLyC, 
from stellar mass to total number of LyC photons produced, is regulated by the IMF and 
model atmospheres. The user has the option of controlling these parameters to determine the 
resulting critical SFR density, subject to several observational constraints. This simulator 



can be accessed at http://casa.colorado.edu/~harnessa/SFRcalculator with an html 
interface for easy use. 

Two standard tests of the SFR use the simulator to calculate the ionization histories 
of H II and He III and compare them to the ionized volume filling factor, Qmi{z), and the 
CMB optical depth, Te{z). The calculator derives Te{z) by integrating the differential form 
of Eq. (3) for various SFR histories and parameters. The average evolution of Qmi is found 
by numerical integration of the rate equation (Madau, Haardt, & Rees 1999) expressing the 
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sources and sinks of ionized zones, 

dQmi _ ni.yc 
dt (uh) 

The source term, riLyC = Psfr/cscQlyC, represents the net ionizing photon production rate, 
computed from selected models of the SFR density, /?sfr- Here, {uh) is the mean hydrogen 
number density, and (tree) is the hydrogen recombination timescale, which depends on Ch 
as shown in Eq. (6). The IGM is assumed to be fully ionized when Qhii = 1- We integrate a 
similar equation for QhoIii to follow He II photoionization in the QSO 4-ryd continuum. We 
use the QSO emissivities at 1 ryd (Haardt & Madau 2012), extrapolated to 4 ryd assuming 
a spectrum with specific flux F^, oc u~^'^. 

Figure 5 shows ionization histories, quantified by Qhii(^) and Te{z), for various values of 
clumping factor, escape fraction, IGM temperature, and LyC-production efficiencies (Qhyc)- 
We adopt the SFR history from Trenti et al. (2010) with an evolving luminosity function. Our 
standard model adopts Ch = 3, /esc = 0.2, T4 = 2, and QljC = 0.004 (4x 10^° photons/M©). 
Different curves show the effects of changing Ch and /esc, including two models in which 
these parameters evolve with redshift. With these parameters, we are typically able to 
complete reionization by z ^ 7, consistent with observations of Gunn-Peterson troughs, 
redshift evolution in Lya emitters, and IGM neutral fraction. The model with a constant 
/esc = 0.05 does not complete the ionization until z ^ 4, which is far too late. Thus, we be 
believe that LyC escape fractions must be considerably larger (> 20%) perhaps evolving to 
higher values at 2; > 6 shown by green and magenta curves. 

Figure 6 compares the dependence of Qmiiz) and Te{z) on SFR histories, computed 
with an evolving luminosity function (Trenti et al. 2010) and new calculations (Haardt & 
Madau 2012) of the star formation rate density, Psfr{z). Figure 7 compares the effects of 
two choices of model atmospheres, with a fixed SFR history from Trenti et al. (2010). The 
on-line calculator provides ionization fractions, Qmi{z) and Qneiii{z), together with CMB 
optical depth, Te{z). In the simulator, users can select other IGM parameters and SFR 
histories. The main difference between the two SFR models lies in the assumptions at z > 8. 
Haardt & Madau (2012) rely on an empirical extrapolation of the star formation rate as a 
function of redshift, while Trenti et al. (2010) adopt a physically motivated model based on 
the evolution of the dark-matter halo mass function. The two approaches are similar at 2; < 8, 
but differ significantly at higher redshift, where an empirical extrapolation does not capture 
the sharp drop in the number density of galaxies observed at 2; ~ 10 (see Figure 8 in Oesch 
et al. 2011). As a consequence, the reionization history from the Haardt & Madau (2012) 
model is more extended at high z, especially when the efficiency of reionization is increased 
because of evolving clumping factor and escape fraction (our preferred models, green and 
magenta curves in Fig. 5). The two models yield quite different predictions for the duration 
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of reionization, defined as tlie redsliift interval, Az, over whicli Qmi evolves from 20% to 80% 
ionized. Our preferred SFR models (green and magenta lines in Fig. 5) have Az ~ 3 (from 
z ~ 10.5 to z ~ 7.5), whereas the Haardt-Madau SFR histories exhibit a more extended 
interval, Az ~ 6 (from z ^ 13 to z ^ 7). This difference in Az could be tested by upcoming 
21-cm experiments; see Bowman & Rodgers (2010) for an initial constraint, Az > 0.06. 

Figure 8 illustrates a third constraint on the reionization epoch (Pritchard et al. 2010; 
Lidz et al. 2011) comparing SFR histories with estimates of the ionizing background at 
^ = 5.5 ± 0.5 (Fan et al. 2006; Bolton & Haehnelt 2007; Haardt & Madau 2012). The 
LyC co-moving emissivity (in photons s~^ Mpc~^), defined as nLyC = PsFR/escQLyC, can be 
related to the ionizing background at 2; = 5 — 6, using recent estimates of the hydrogen 
photoionization rate, r}ii[z), from Haardt & Madau (2012) and the LyC mean free path, 
Ahi, from Songaila & Cowie (2010). The LyC emissivity is proportional to the star formation 
rate density, computed from our halo mass function model (Trenti et al. 2010), integrated 
down to absolute magnitudes Mab = —18 (Bouwens et al. 2011a) or to Mab = —10, the 
faint limit suggested by Trenti et al. (2010). We adopt a fiducial LyC production parameter 
QhyC = 0.004, corresponding to 4 x 10^*^ LyC photons produced per Mq of star formation, 
and we use two different models for LyC escape fraction, /esc (constant at 20% and varying 
with redshift). For quantitative values, we assume an ionizing background with specific 
intensity, = Jo{v/vo)~'^ (in units erg cm~^ s~^ sr~^ Hz~^) with power-law index a 2 at 
energies above huQ = 1 ryd. The hydrogen photoionization rate is Fhi = [4:njQao/h{a + 3)] 
for a hydrogen photoionization cross section a,^ ^ aQ^u/uQ)'^ with ctq = 6.3 x 10^^^ cm^. 
For this spectrum, the frequency-integrated ionizing intensity is Jtot = </oi^o/('^ — 1), and the 
LyC photon flux (photons cm~^ s~^) integrated over all solid angles is $LyC = {4:7i Jq / ha) . 
Because we analyze the photon mean-free path, we normalize Jq to $LyC and Fhi, 

/za$LyC _ h{a 3) „ 



Air 



HI 



By approximating $LyC as the product of the LyC emissivity and mean-free path, we arrive 
at the calibrations: 
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Finally, we relate $LyC to the star-formation rate Psfr and ionization rate Fhi 
PSFR = (0.056 Mq yr-i Mpc 
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Here, we have scaled the parameters to the same values assumed by Lidz et al. (2011), namely 
Fhi = 5 X 10~^^ s~^, /esc = 0.1, an ionizing spectral index a = 2, and redshift z = 5.5, at 
which Songaila & Cowie (2010) fit Ahi ~ 9.8 proper Mpc. Our parameter QljC = 0.004 
corresponds to the Lidz et al. (2011) LyC photon production calibration, 10^^'^ s~^ per 
Mq yr~^ of star formation. However, our coefficient, 0.056 Mq yr~^ Mpc~'^, is slightly larger 
than their value, 0.039 Mq yr~^ Mpc~^, an effect that may arise from our different method 
of relating Fhi to the ionizing radiation field and SFR. 

More careful examination of these parameters suggests that, at 2; ~ 5.5, the hydrogen 
ionization rate Fhi ~ 3.6 x 10~^^ s~^, given in Table 3 of Haardt & Madau (2012), and the 
observed mean free path, Ahi, at z = 5.5 may be closer to 6 proper Mpc (see Figure 10 of 
Songaila & Cowie 2010). Rescaling to those two parameters, we find a similar coefficient of 
0.066 Mq yr-i Mpc"l This SFR density at z = 5.5 ± 0.5 (see also Figure 8) is comparable 
to the critical value needed to maintain reionization at 2 = 7, but the observed rates and 
mean free paths are declining rapidly with redshift. All three constraints on SFR suggest 
that full reionization is more likely to occur at redshift Zj-d ~ 7 than at 2rci = 10. 



4. DISCUSSION AND IMPLICATIONS 

The major results of our study can be summarized with the following points: 

1. We calculated the critical star formation rate required to maintain a photoionized 
IGM, incorporating four free parameters {Ch, /esc, Tf>, Qljc) that control the rates of 
LyC photon production and radiative recombination. Our best estimate at 2; = 7 is 
Pcrit ~ (0.018 Mq yr-1 Mpc-3)[(1 + z)/8f (Ch/S) {0.2/ U,)T^°-^^^ for fiducial values of 
IGM clumping factor Ch = 3, LyC escape fraction /esc = 0.2, temperature Tg = 10"^ K, 
and standard IMFs and low-met allicity stellar atmosphere {QljC = 0.004). 

2. An epoch of full reionization at Zj-^ ~ 7 is consistent with recent optical/IR measure- 
ments of SFR history and a rising IGM neutral fraction at 2 = 6 — 8, marking the 
tail end of reionization. These observations include the decrease in numbers of high- 2; 
galaxies and Lya emitters and IGM damping-wing intrusion into the Lya transmission 
profiles in high- 2; QSOs. 

3. Our newly developed SFR and reionization calculator, now available on-line at 



http : //casa. Colorado . edu/~harnessa/SFRcalculator , allows users to calculate the 



ionization history, critical SFR, and CMB optical depth, and to assess whether observed 
SFRs, IMFs, and other parameters are consistent with IGM reionization. 
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4. Reconciling late reionization at z^^i ~ 7 with Tg = 0.088 ± 0.015 of the CMB likely 
requires an epoch of partial ionization. A fully ionized IGM back to z^^i = 7 produces 
Tg ~ 0.050, and an additional optical depth, Arg ~ 0.01 — 0.04, may arise from early 
sources of UV/X-ray photons at z > 7. Alternatively, one can appeal to the likelihood 
contours from WMAP, which allow optical depths as low as 0.06 (95% CL.) The Planck 
experiment may clarify the situation in several years. 

5. If the EoR is more complex, as we suggest, then redshifted 21-cm experiments should 
focus on the interval 7.5 < z < 10.5, corresponding to frequencies 124-167 MHz, for the 
maximum signal of IGM heating (T;, 27 mK) produced when the hydrogen neutral 
fraction xhi ~ 0.5 (Pritchard et al. 2010). 

We conclude by speculating about which future observations can best constrain the EoR. 
Additional data from from the Planck mission^ should provide confirmation of the CMB 
optical depth with smaller error bars. This information will constrain the additional amount 
of ionization at 2; > 7 and narrow the range Arg = 0.01 — 0.04 produced by high-redshift 
sources in the partially ionized IGM. Ongoing surveys for high-z galaxies, Lya emitters, and 
QSO near-zone sizes will better quantify the rise of neutral fraction, xhi at 2; > 6.5. On the 
theoretical front, we are running larger IGM simulations on 1536^ and 2048^ unigrids and 
will add discrete sources of ionizing radiation and radiative transfer in order to capture the 
heating and clumping more accurately. We will also include source turn-on at 2 > 9 to test 
the decrease and recovery of Ch, as seen in Figures 3 and 4. Finally, we plan to carry out more 
detailed modeling of the H I (21-cm) signal, coupled to the kinetic and spin temperatures 
driven by heating at 2; > 7. As described earlier, the duration of the reionization transition 
and the 21-cm emission during the heating phase could provide discriminants of various SFR 
histories at 2; > 6. 

This work was supported by grants to the Astrophysical Theory Program (NNX07- 
AG77G from NASA and AST07-07474 from NSF) at the University of Colorado Boulder. We 
thank Joanna Dunkley, Chris Carilli, Richard Ellis, and Piero Madau for useful discussions on 
reionization and CMB optical depth statistics. We are grateful to the referee for suggesting 
additional model comparisons with the ionizing background at 2 = 5 — 6. 



^See mission website |http : //www . rssd . esa . int /Planck 
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Table 1. Sample values of QljC and pcrit for different IMFs 



M^in (Mq) 


A/max (Mq) 


a 


Z 






Pent 


Pent 


0.1 


100 


2.35 


Zq 


0.00236 


0.00286 


0.0306 


0.0253 


0.1 


100 


2.35 


O.O2Z0 


0.00397 


0.00558 


0.0181 


0.0129 


0.1 


100 


2.35 





0.00401 


0.00752 


0.0180 


0.0096 


0.1 


100 


2.35 


0.2Zq 


0.00365 




0.0197 




0.1 


100 


2.35 


O.IZq 


0.00383 




0.0188 




1.0 


100 


2.35 


O.IZq 


0.00976 




0.0074 




0.1 


100 


2.00 


O.IZq 


0.01267 




0.0057 




0.1 


200 


2.35 


O.IZq 


0.00543 




0.0133 





''Sutherland & ShuU unpublished model atmospheres 
'^Schaerer model atmospheres 



Note. — Production efficiency of ionizing (LyC) radiation, QhyC, in units of lO"'' 
photons/MQ. Critical SFR, Pcriti is given in units of Mq yr^-"^ Mpc^'^ (comoving) as- 
suming Ch = 3, /esc = 0.2, and T4 = 1. 



Table 2. Parameters of Simulations 



Run 


I 

(h-i Mpc) 


7V^/f 

cells 


^uv 


Relative Strength 
of Background 


50.1024_7'' 


50 


1024 


7 


1 


50.512.7'^ 


50 


512 


7 


1 


50_512_9 


50 


512 


9 


1 


50_512_9.2 


50 


512 


9 


2 


50.512.0 


50 


512 


N/A 





50.1536.9 


50 


1536 


9 


1 



''•''Simulations 50.1024.2, 50.512.2 from Smith et al. (2011). 

Note. — zyy is the redshift at which the ionizing background 
radiation (Haardt & Madau 2001) is turned on, ramping up to a 
constant value hy z = zyv — 1 
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Fig. 1. — Comparison of two methods for calculating the clumping factor for simulation 
50_1024_7 (UV background turned on at 2 = 7) with different weights. Red lines: clumping 
calculated via density field (DF) method. Blue lines: clumping via recombination rate (RR) 
method (see Section 2.2). Solid lines correspond to weighting by volume, and dashed lines 
correspond to weighting by baryon overdensity. We believe the RR method, with volume 
weighting, is a more accurate measure of clumping. 
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Fig. 2. — Effects of different "density cuts" on the clumping factor from our 1536'^ simulation, 
using the RR method in Equation (15). Solid line: clumping factor summed over cells with 
overdensity between 1 < A5 < 100, excluding low-density voids with Af, < 1. Dashed line: 
including all cells with A(, < 100. 
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Fig. 3. — Comparison of convergence and cosmic variance among simulations. Values of Ch 
agree for 1024'^ and 1536'^ simulations, with variance computed over eight sub-volumes in 
our simulations. Top panel: eight 768'^ sub-volumes in 1536^ run. Bottom panel: eight 512'^ 
sub-volumes in 1024^ run. In each panel, our original 512'^ run is shown as dotted black line, 
and values for 1536^ and 1024"^ runs as heavy black lines. Large variations in Ch high 
redshifts arise from small numbers of cells prior to turn-on of the ionizing background at 
z = 9 (top) and z = 7 (bottom). For the 1536^ simulation (top panel), the clumping factor 
is well-fitted by Cniz) = (2.9) [(1 + z)/6]^^-^ between 5 < 2 < 9. 
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Fig. 4. — Clumping factor calculated via the density field (DF) method computed from three 
512^ IGM simulations, with sums weighted by volume (see eq. [14]). Red solid line has UV 
background (Haardt & Madau 2001) turned on at 2; = 9 (simulation 50_512_9), blue dashed 
line at 2; = 7 (simulation 50_512_7), and green dot-dashed line has no background (simulation 
50_512_0). We observe two tracks for Ch{z): a high-track for redshifts above turn-on of the 
UV background, and a lower track after turn-on. During most of the reionization epoch, 
from 6 < z < 12, Cr lies between 1.5 and 3. 
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Fig. 5. — History of ionization fraction Qmii^) and CMB optical depth Te{z) versus redshift 
(eqs. [3] and [17]) computed with Te = 2 x 10^ K, QLyC = 0.004 (4 x 10^° photons per M©) 
and using SFR history and evolving luminosity function from Trenti et al. (2010) integrated 
down to A^AB = —10. The IGM is assumed to be fully ionized when Qmi = 1- Blue line: 
constant Ch = 3 and /esc = 0.05 (this model reionizes too late to be viable). Red line: 
constant Ch = 3 and /esc = 0.2. Green line: constant Ch = 3 but redshift-dependent 
/esc = 1-8 X 10-'^(1 + z)3-^ (Haardt & Madau 2012). Magenta line: redshift-dependent 
Ch = 1 + 432-1-^1 (Pawlik et al. 2009) and Uc = 1-8 x 10-^(1 + zf-^. Our two preferred 
models with variable Ch or /esc (green, magenta) naturally produce Zj-d ~ 7. Solid black 
circles indicate redshifts of 20%, 50%, and 80% ionization; the duration Az ^ 3 is defined 
between 20% and 80% points. 
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Fig. 6. — Top: Volume filling factors, Qmi, with same color codes as Fig. 5, comparing two 
models for SFR history. Solid lines show SFR history from Trenti et al. (2010), and dashed 
lines show SFR from equation (53) in Haardt & Madau (2012), which is larger and more 
extended at higher redshifts, particularly for z > 8. Bottom: Corresponding CMB optical 
depths Te{z). Additional optical depth may arise from sources at z > 7. 
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Fig. 7. — Example of the output for Qmiiz), QHeiii(^), and Te{z) from our on-line 
SFR/reionization simulator. Solid lines show model atmospheres labeled as Sutherland & 
ShuU (S&S), while dashed lines are for Schaerer (2002, 2003), assuming Z = 0.02Zq. Mod- 
els assume SFR history of Trenti et al. (2010) and parameters with constant Ch = 3 and 
redshift-dependent /esc = 1-8 x 10~^(1 + zf-'^ (Haardt & Madau 2012). Dot-dashed line 
shows He III ionization history, Qucuii^), computed for QSO (4-ryd continuum) ionization 
as described in Section 3.3. In these models, H I is fully ionized by 2; ~ 7 and He III by 
z ~ 2.7. Additional optical depth may arise from sources at 2; > 7. 
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Fig. 8. — Top: Co-moving SFR density (psfr) from halo mass function model (Trenti et al. 
2010) integrated for galaxies down to Mab = —18 (solid blue line) and to —10 (dotted 
green line). Red data points are from Bouwens et al. (2011a) for galaxies down to Mab = 
— 18. Solid data point at 2 ~ 5.5 is from ionizing background constraint (Eq. 21), using 
Haaradt & Madau (2012) model of Th{z) and LyC mean free paths from (Songaila & Cowie 
2010). Bottom: Co- moving emissivity, nLyC? of LyC photons, corresponding to above SFRs, 
integrated to Mab = — 18 and —10. We use a LyC production parameter QljC = 0.004 
^2q53.i photons s~^ per Mq yr~^) and two different escape-fraction models: /esc = 0.2 
(blue curves and data points) and fcsc{z) = (1.8 x 10~^)(l-|-2;)^'^ (red curves and data points). 



